From LoCoH to T-LoCoH
Eric Dougherty & Dana Seidel
November 29, 2017
Our walkthrough today will generally follow the guideline laid out by Andy Lyons (an alum of the Getz lab) when he first published his T-LoCoH paper in Movement Ecology. The PDF of that tutorial can be found in our Day 4 folder and also at: T-LoCoH Tutorial
Just like we have been doing for the past few days, the first step will be installing the tlocoh package and getting it loaded into our R session. Unlike some of our previous examples, though, this particular package is housed on R-Forge rather than CRAN. This means we will need to use an additional argument or two:
#install.packages("tlocoh", dependencies=T, repos=c("http://R-Forge.R-project.org"))
library(tlocoh)## Loading required package: sp## Warning: package 'sp' was built under R version 3.3.2## T-LoCoH for R (version 1.40.04)
## URL: http://tlocoh.r-forge.r-project.org/
## Bug reports: tlocoh@gmail.comThe first step for us will be to load in one of the datasets included in the package. All we need to do is:
data(toni)Toni is the name that researchers from the Getz Lab gave to one of the African buffalo they studied in the Kruger National Park in South Africa during a study of bovine tuberculosis. Let’s take a look at her data using the base::head command:
dim(toni)## [1] 6371 4head(toni)## id long lat timestamp.utc
## 17930 toni 31.75345 -24.16950 2005-08-23 06:35:00.000
## 17931 toni 31.73884 -24.15402 2005-08-23 07:34:00.000
## 17932 toni 31.73969 -24.15359 2005-08-23 08:34:00.000
## 17933 toni 31.73874 -24.15329 2005-08-23 09:35:00.000
## 17934 toni 31.73946 -24.15336 2005-08-23 10:34:00.000
## 17935 toni 31.73898 -24.15363 2005-08-23 11:35:00.000We can see that we have 6371 data points for Toni. We can also see that the positional fixes are in the form of latitude and longitude. In order to use these data for our analyses here, we will need to convert these lat-long points into a meaningful geographic projection such as UTM, which will transform the points from degrees to meters. Knowing that these data are from South Africa, we can determine the UTM zone in which these data were collection (Zone 36S). In order to reproject the data, we will need two more packages that we have dealt with in the past, sp and rgdal. The first order of business is to create a SpatialPoints object with a defined projection (longlat, in this case). Then, we will transform the underlying projection to one more suited to our needs:
library(sp)
library(rgdal)## Warning: package 'rgdal' was built under R version 3.3.2## rgdal: version: 1.2-8, (SVN revision 663)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
## Linking to sp version: 1.2-5toni.sp.latlong <- SpatialPoints(toni[ , c("long","lat")], proj4string=CRS("+proj=longlat +ellps=WGS84"))
toni.sp.utm <- spTransform(toni.sp.latlong, CRS("+proj=utm +south +zone=36 +ellps=WGS84")) Now we can take a look at the coordinates to make sure they are no longer in degrees:
toni.mat.utm <- coordinates(toni.sp.utm)
head(toni.mat.utm)## long lat
## 17930 373372.0 7326443
## 17931 371872.3 7328144
## 17932 371958.1 7328192
## 17933 371861.3 7328225
## 17934 371934.5 7328218
## 17935 371886.0 7328187Looks like we were successful! The only thing left to do is change the column labels from ‘long’ and ‘lat’ to ‘x’ and ‘y’
colnames(toni.mat.utm) <- c("x", "y")Now we have a matrix (with two columns) of all of the positional fixes obtained from Toni during the study. One more thing we will need to check is the timestamp on these points. We have already seen the POSIX class in action, but one of the important arguments for defining that class is the timezone in which the points were taken. While some researchers collect their data in the local timezone, others prefer to standardize to UTC. In this case, the latter was selected, and we can tell by the name of the column in the original toni data file! Let’s use the base::as.POSIXct command to define the time zone and the class of the timestamps:
toni.gmt <- as.POSIXct(toni$timestamp.utc, tz="UTC")Now that the timezone has been set, we may want to transform it into the time zone where the animal was actually located, as the time of day (as experienced by the animal) may be important for our understanding of their movement. For example, it would be really difficult to determine if an animal is nocturnal or diurnal or crespuscular (active at sunrise and sunset) if you don’t know when sunrise and sunset are relative to the time at which the points were collected. In this case, POSIX recognizes the timezone of Kruger as “Africa/Johannesburg”
local.tz <- "Africa/Johannesburg"
toni.localtime <- as.POSIXct(format(toni.gmt, tz=local.tz), tz=local.tz)We now have everything we need (in the correct formats) to move on to an actual LoCoH analysis!
The first step in using LoCoH is the creation of the fundamental LoCoH object: lxy. To do this, we use the tlocoh::xyt.lxy command. This unique object requires a few inputs, including the UTM coordinates (and the associated projection), the timestamps of the points, and the ID (name) of the individual.
toni.lxy <- xyt.lxy(xy=toni.mat.utm, dt=toni.localtime, id="toni", proj4string=CRS("+proj=utm +south +zone=36 +ellps=WGS84"))## 595 duplicate xy-time-id rows removedIt turns out there were 595 instances in the dataset of 6371 points where the timestamp and the coordinates were duplicated. In the creation of the lxy object, these points were removed, leaving a total of 5776 points. You could tell this from digging into our new lxy object which currently consists of 5 elements (in the form of a list): pts, dt.int, rw.params, anv, and comment.
nrow(toni.lxy$pts)## [1] 5776nrow(toni.lxy[[1]])## [1] 5776The structure of the lxy object allows for various aspects of it to be inspected in several different ways. These two commands told us the same information, but one accessed the first element of the list directly by name whereas the other called it without calling for the pts element by name.
We can also take a look at some of the other descriptors of the object:
summary(toni.lxy)## Summary of LoCoH-xy object: toni.lxy
## ***Locations
## id num.pts dups
## toni 5776 9
## ***Time span
## id begin end period
## toni 2005-08-23 2006-04-23 243.3 days
## ***Spatial extent
## x: 369305.5 - 391823.9
## y: 7305737.9 - 7330491.3
## proj: +proj=utm +south +zone=36 +ellps=WGS84
## ***Movement properties
## id time.step.median d.bar vmax
## toni 3600 (1hs) 173.7575 0.9267969
## ***Ancilliary Variables:
## -none-
## ***Nearest-neighbor set(s):
## none savedThis illustrates the number of points, the length of the movement trajectory, the spatial extent of the animal’s movement, some general properties of the track, and would also have information of ancillary variables and nearest neighbor sets (if these exist, which they do not here). Another way to examine our data would be to tlocoh::plot it using the function in the locoh package.
plot(toni.lxy)
This unique plot function colors the points by the timestamp so that you can get an idea of the recursion in the path (i.e., the way the animal covers the same space at different times). This really helps demonstrate the need for T-LoCoH, but we will get to that in a bit more detail soon. One other feature that we can determine from this plot, which is much more difficult to ascertain when staring at a list of numbers, is that there do not appear to be any significant outliers. This suggests that there are no erroneous points, such as those obtained before the collar was deployed.
Yet another summary of the data can be seen using the tlocoh::hist command, which will display the distribution of locations by date, step length, and sampling interval.
hist(toni.lxy)
The left two panels (Num Locations Over Time and Time Interval) illustrate that the sampling throughout the period was relatively uniform; it does not look like there were too many missed points during the sampling (though there were a few, particualry towards the end). One normally looks for these gaps because we know that missing data is a problem, but in the case of movement trajectories, the opposite (too many points in a short period) may also be a problem. These are called ‘bursts’, and the tlocoh package uses the tlocoh::lxy.plot.freq command to determine whether there are points that were obtained at too high a frequency relative to the expected sampling frequency.
lxy.plot.freq(toni.lxy, cp=T)
This dataset appears to have one little burst of points (the little dot in the lower left corner), so we’ll thin it out. Setting a threshold of 0.2 (meaning that any group of points that are less than 0.2 times the median sampling interval) will be considered a cluster and thinned down to one location.
toni.lxy <- lxy.thin.bursts(toni.lxy, thresh=0.2)Now comes the good part. Presumably, we would like to know something about the home ranging behavior of the animal for which we have data. As we have discussed, there are quite a few ways to do this, but we’ll begin with the LoCoH conception here. In order to understand how we will build the LoCoH home ranges, though, we will need to understand the s parameter first. This parameter is most important when we are creating a T-LoCoH based home range, but for our purposes here, we will still need to set the parameter, so we will dicsuss it briefly. The magnitude of s indicates the degree to which local hulls are local in time as well as space. Essentially, it serves to transform distance into time-scaled distance, where the temporal separation between points is considered in addition to the spatial separation. As I mentioned, it is not especially important just yet becasue we are going to set it at zero for our first example, but it will become a key component in later examples.
So let’s begin by constructing a model of space-use without considering time. The first thing we will want to do is identify nearest neighbors. The nearest neighbors can be defined in several ways, but the simplest is the k method, which just searches around each point for the k nearest points. Because we are not considering time yet (s=0), these nearest points will be those that are closest in space. Let’s start out by selecting k=25. We will add these to our existing lxy object using the command tlocoh::lxy.nn.add:
toni.lxy <- lxy.nn.add(toni.lxy, s=0, k=25)## Finding nearest neighbors for id=toni (n=5775),
## num.parent.pts=5775, mode=Fixed-k, k=25, s=0, method=Euclidean
## - computing values of kmax, rmax, and amax...Done
## - set of neighbors (re)named: toni|vmax|s0|n5775|kmax25|rmax57.8|amax701
##
## Done. Nearest neighbor set(s) created / updated:
## toni|vmax|s0|n5775|kmax25|rmax57.8|amax701
## Total time: 0.5 secsNow, if we take a look at our summary of the lxy object, you will see this new nearest neighbor set listed:
summary(toni.lxy)## Summary of LoCoH-xy object: toni.lxy
## ***Locations
## id num.pts dups
## toni 5775 9
## ***Time span
## id begin end period
## toni 2005-08-23 2006-04-23 243.3 days
## ***Spatial extent
## x: 369305.5 - 391823.9
## y: 7305737.9 - 7330491.3
## proj: +proj=utm +south +zone=36 +ellps=WGS84
## ***Movement properties
## id time.step.median d.bar vmax
## toni 3600 (1hs) 173.7452 0.9267969
## ***Ancilliary Variables:
## -none-
## ***Nearest-neighbor set(s):
## 1 toni|vmax|s0|n5775|kmax25|rmax57.8|amax701As we would expect, the kmax listed there is 25 (we set it at 25, so it better be 25!), but you’ll also see an rmax and an amax value. These can be used for alternative search methods for identifying the nearest neighbors. In the r-method, you set a radius around each point, and all of the points within that radius are treated as the nearest neighbors. Unlike the k-method, this will result in a nearest neighbor set where each point has a different number of nearest neighbors. The a-method is similar in that it is likely to result in different numbers of nearest neighbors for each point, but for this method, a total distance is defined (say, 701) and the algorithm selects the closest neighboring point and then subtracts that distance (say 50) from the total (now, 651). The algorithm continues adding neighbors and subtracting their distances to the point in question from the total distance. When the next closest neighbor would require more distance than remains in the total distance, it is ignored and the other neighbors represent that point’s nearest neighbor set. The maximum r and a values indicate that this nearest neighbor set defined using the k-method would be sufficient for r or a values up to those maxima without rerunning the above command.
Now that we have expended all of this energy making this lxy object, we might want to save it to our working directory using the tlocoh::lxy.save command:
#lxy.save(toni.lxy, dir=".")The building blocks of all T-LoCoH analyses are hulls, which are simply minimum convex polygons constructed around each point from a set of nearest neighbors. Since we’ve already identified 25 nearest neighbors for each point, we can create hulls with up to 25 nearest neighbors each. We will use the tlocoh::lxy.lhs command to create an lhs object, the next in our progression.
toni.lhs <- lxy.lhs(toni.lxy, k=3*3:8, s=0)## Using nearest-neighbor selection mode: Fixed-k
## Constructing hulls and hull metrics...
## toni: 9 duplicate points were randomly displaced by 1 map unit(s)
##
## toni.pts5775.k9.s0.kmin0
## Found a suitable set of nearest neighbors
## Identifying the boundary points for each parent point
## Converting boundary points into polygons
## Calculating area and perimeter...Done.
## Calculating the time span of each hull...Done.
##
## toni.pts5775.k12.s0.kmin0
## Found a suitable set of nearest neighbors
## Identifying the boundary points for each parent point
## Converting boundary points into polygons
## Calculating area and perimeter...Done.
## Calculating the time span of each hull...Done.
##
## toni.pts5775.k15.s0.kmin0
## Found a suitable set of nearest neighbors
## Identifying the boundary points for each parent point
## Converting boundary points into polygons
## Calculating area and perimeter...Done.
## Calculating the time span of each hull...Done.
##
## toni.pts5775.k18.s0.kmin0
## Found a suitable set of nearest neighbors
## Identifying the boundary points for each parent point
## Converting boundary points into polygons
## Calculating area and perimeter...Done.
## Calculating the time span of each hull...Done.
##
## toni.pts5775.k21.s0.kmin0
## Found a suitable set of nearest neighbors
## Identifying the boundary points for each parent point
## Converting boundary points into polygons
## Calculating area and perimeter...Done.
## Calculating the time span of each hull...Done.
##
## toni.pts5775.k24.s0.kmin0
## Found a suitable set of nearest neighbors
## Identifying the boundary points for each parent point
## Converting boundary points into polygons
## Calculating area and perimeter...Done.
## Calculating the time span of each hull...Done.
## The following hullsets were generated:
## toni.pts5775.k9.s0.kmin0
## toni.pts5775.k12.s0.kmin0
## toni.pts5775.k15.s0.kmin0
## toni.pts5775.k18.s0.kmin0
## toni.pts5775.k21.s0.kmin0
## toni.pts5775.k24.s0.kmin0
## Total time: 19.9 secsWhat we are doing here is creating 6 different hullsets with k values of 9, 12, 15, 18, 21, and 24 (note we are still setting s=0 because we are not considering the temporal component yet). It is useful to create a selection of hullsets so that we can compare them and choose the optimal one for our purposes. The downside of creating all of these extra hullsets is that it may take a little while to run. Once it is finished, we can get an idea of what areas are included in our hullsets using the tlocoh::plot command again. This time, we will specify that we want to plot the hulls and we’ll add all of the plots to a single page for easier comparison.
plot(toni.lhs, hulls=TRUE, figs.per.page=6)
That looks pretty good, but it doesn’t tell us too much just yet. Let’s create isopleths for our hullset. Isopleths are aggregations of hulls sorted in such a way as to reveal something about space use. The default settings for tlocoh::lhs.iso.add sorts hulls according to density, so the isopleths reflect the likelihood of occurrence, which is a proxy for intensity of use. Then we can plot these (using the iso=TRUE specification instead of hulls=TRUE) to see how they compare to one another and to the hullset plots.
toni.lhs <- lhs.iso.add(toni.lhs)## Merging hulls into isopleths
## toni.pts5775.k9.s0.kmin0
## Sorting hulls by area...Done.
## Unioning hulls
##
|
| | 0%
|
|++++++++++ | 20%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++| 100%
## 5 invalid polygon(s) removed
## toni.pts5775.k12.s0.kmin0
## Sorting hulls by area...Done.
## Unioning hulls
##
|
| | 0%
|
|++++++++++ | 20%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++| 100%
## 8 invalid polygon(s) removed
## toni.pts5775.k15.s0.kmin0
## Sorting hulls by area...Done.
## Unioning hulls
##
|
| | 0%
|
|++++++++++ | 20%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++| 100%
## 2 invalid polygon(s) removed
## toni.pts5775.k18.s0.kmin0
## Sorting hulls by area...Done.
## Unioning hulls
##
|
| | 0%
|
|++++++++++ | 20%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++| 100%
## toni.pts5775.k21.s0.kmin0
## Sorting hulls by area...Done.
## Unioning hulls
##
|
| | 0%
|
|++++++++++ | 20%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++| 100%
## toni.pts5775.k24.s0.kmin0
## Sorting hulls by area...Done.
## Unioning hulls
##
|
| | 0%
|
|++++++++++ | 20%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++| 100%
## 1 invalid polygon(s) removed
## Total time: 10.3 secsplot(toni.lhs, iso=TRUE, figs.per.page=6)
According to the tutorial, we are looking to select the k value that leads to a set of hulls where the heavily used area doesn’t look like Swiss cheese, but also doesn’t cut across unused areas. This is essentially finding a balance between Type I error (including area that isn’t part of the home range) and Type II error (omitting area the animal used). This is a somewhat subjective choice and it will depend on the particular research question. Once we select one that we think satisfies those criteria, we can take a closer look at our choice before moving forward with our analyses. Using the allpts=TRUE argument, we can take a look at what the points are included in our selected home range.
plot(toni.lhs, iso=T, k=15, allpts=T, cex.allpts=0.1, col.allpts="gray30", ufipt=F)
In addition to the isopleth plots, we can use two other metrics to verify that our choice of k value is sufficient. The first is the isopleth area curves, which will plot the area included in a set of different isopleth levels (ranging from 0.15 to 0.95 by 0.10). We want to check these plots to make sure that there are no sharp jumps between k values that indicate that a relatively small change in k results in a large increase in included area (likely a false commission). To look at this plot, we use the tlocoh::lhs.plot.isoarea command:
lhs.plot.isoarea(toni.lhs)
There don’t appear to be any major jump (i.e., the curves are all resonably smooth, no matter the k value). This means we shouldn’t rule out any of the k values we have tested. The next metric we could look at is the edge to area ratio. A simple plot of this measure across the same set of isopleths can be created using the tlocoh::lhs.plot.isoear command. We are looking to exclude k values that result in very high edge:area ratios, which are indicative of the Swiss cheese pattern (i.e., many small holes). This is particularly important at some of the isopleths associated with the core area (0.35, 0.45, 0.55) where we would expect relatively few holes.
lhs.plot.isoear(toni.lhs)
We can see that a k of 9 results in some pretty Swiss-cheesy looking isopleths at the lower end of the spectrum, but by k=15, these are largely filled in. So, if k=15 still looks sufficient based on these metrics and the isopleth map, we can use the tlocoh::lhs.select command to select this value for subsequent analyses:
toni.lhs.k15 <- lhs.select(toni.lhs, k=15)Now that we have this lovely hullset, let’s make sure we save it so that we don’t need to rerun all of these steps:
lhs.save(toni.lhs.k15)## LoCoH-hullset toni.lhs.k15 saved as 'toni.lhs.k15' to:
## /Users/ericdougherty/Desktop/GitHub/MovEco-R-Workshop/Materials/Day4/toni.n5775.s0.k15.iso.lhs.03.RDataLet’s take a few steps back and think about incorporating the temporal component to our analysis. As we saw in that first plot of our track, the animal definitely used the same space at different times. The current lhs object that we created ignored the separation of these points in time and creates a hull out of the nearest points in terms of spatial proximity alone. But knowing what we do, we should really consider time and at least examine how things change.
We are going to return to our lxy object (which already has one set of nearest neighbors attached), and create a new nearest neighbor set with a non-zero s parameter. But we don’t want to select this value willy-nilly! Fortunately, Andy has a few different methods that we can use to obtain a reasonable s value. The first is to pick an s value such that 40-80% of the hulls are ‘time-selected.’ All you need to do for this method is use the tlocoh::lxy.ptsh.add command:
toni.lxy <- lxy.ptsh.add(toni.lxy)## id: toni
## Randomly selected 200 of 5775 points
## Finding 10 nearest neighbors for 200 sample points
## Finding s for ptsh=0.98
## s=0.005, s=0.01, s=0.02, s=0.04, s=0.08, s=0.16,
## Finding s for ptsh=0.1 (+/- 0.01)
## s=0.0025, s=0.00125, s=0.000625,
## Finding s for ptsh=0.2 (+/- 0.01)
## s=0.001875,
## Finding s for ptsh=0.3 (+/- 0.01)
## s=0.0075, s=0.00625, s=0.005625, s=0.0053125,
## Finding s for ptsh=0.4 (+/- 0.01)
## Finding s for ptsh=0.5 (+/- 0.01)
## s=0.015, s=0.0125, s=0.01375,
## Finding s for ptsh=0.6 (+/- 0.01)
## s=0.0175,
## Finding s for ptsh=0.7 (+/- 0.01)
## s=0.03, s=0.025,
## Finding s for ptsh=0.8 (+/- 0.01)
## s=0.035,
## Finding s for ptsh=0.9 (+/- 0.01)
## s=0.06, s=0.05,
##
## Done with toni
This creates a set of 200 sample points, each with 10 nearest neighbors. This algorithm goes through a set of different s values and determines what proportion of hulls are time-selected. Though this is a somewhat stochastic process, we can see that s values between approximiately 0.01 and 0.035 fall between proportions of time-selected hulls ranging from about 40-80%. For simplicity, there is also a nice plot that shows the change in the proportion of time selected hulls over various s values to find an optimal point. It seems like 0.03 is a reasonable choice.
The second method is a little more involved. First we need to select a time interval of interest, which will be based on the temporal scale of the question we are asking. For example, we may be interested in daily foraging patterns. If that is the case, we may not want to treat points that occur more than 24 hours apart in time as nearest neighbors of one another, even if they occur close in space. If we are thinking about seasonal patterns, however, points that occur several weeks apart may still reasonably be considered nearest neighbors. If you are unsure exactly what time interval you are most interested in, you can use the tlocoh::lxy.plot.pt2ctr command which plots the distance of each point to the centroid of the entire dataset. This may help isolate the ‘natural’ frequencies in the data.
lxy.plot.pt2ctr(toni.lxy)
Looking at that, it is a little difficult to pick out especially notable patterns, but depending on your data, something may jump out. Even without using that, we can call the `tlocoh::lxy.plot.sfinder command:
lxy.plot.sfinder(toni.lxy)
This shows that approximately half of all pairs of points will be time-selected at the scale of 24 hours when s is somewhere around 0.005 to 0.01. We can get a closer look by defining a specific set of delta.t values:
lxy.plot.sfinder(toni.lxy, delta.t=3600*c(12,24,36,48,54,60))
That’s better! It looks like the results here support the use of 0.03 as our s value. It would not be especially surprising if the two methods did not converge on a single s value, though. They are approaching the selection of s in very different ways. The fact of the matter is that there is not perfect s value. Changes in that parameter will result in different hullsets that tell us different things. Fortunately, we’ve landed on a value of 0.03 using both methods!
Now that we have an s value chosen, let’s return to our lxy object and create a new set of nearest neighbors. We will use the same k value of 25, but we’ll set an s value this time:
toni.lxy <- lxy.nn.add(toni.lxy, s=0.03, k=25)## Finding nearest neighbors for id=toni (n=5775),
## num.parent.pts=5775, mode=Fixed-k, k=25, s=0.03, method=TSD:vmax
## - computing values of kmax, rmax, and amax...Done
## - set of neighbors (re)named: toni|vmax|s0.03|n5775|kmax25|rmax1301.9|amax17017.7
##
## Done. Nearest neighbor set(s) created / updated:
## toni|vmax|s0.03|n5775|kmax25|rmax1301.9|amax17017.7
## Total time: 48.1 secsThis simply adds another nearest neighbor set, which we can see if we use the tlocoh:summary function again:
summary(toni.lxy)## Summary of LoCoH-xy object: toni.lxy
## ***Locations
## id num.pts dups
## toni 5775 9
## ***Time span
## id begin end period
## toni 2005-08-23 2006-04-23 243.3 days
## ***Spatial extent
## x: 369305.5 - 391823.9
## y: 7305737.9 - 7330491.3
## proj: +proj=utm +south +zone=36 +ellps=WGS84
## ***Movement properties
## id time.step.median d.bar vmax
## toni 3600 (1hs) 173.7452 0.9267969
## ***Ancilliary Variables:
## -none-
## ***ptsh s-values computed
## id k n ptsh
## toni 10 200 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
## 0.98
## ***Nearest-neighbor set(s):
## 1 toni|vmax|s0|n5775|kmax25|rmax57.8|amax701
## 2 toni|vmax|s0.03|n5775|kmax25|rmax1301.9|amax17017.7Now there are two items under the nearest neighbor sets! Mission accomplished. We could also have gone about selecting s in a similar fashion as we selected k, essentially creating a set of alternative nearest neighbor sets and comparing the emergent properties of each. In this case, we might want to use the tlocoh::lxy.plot.mtdr command to examine the ratios of diffusion distance to time-scaled distance or tlocoh::lxy.plot.tspan to see the time spans of the resulting hulls.
toni.lxy <- lxy.nn.add(toni.lxy, s=c(0.0003, 0.003, 0.03, 0.3), k=25)## Finding nearest neighbors for id=toni (n=5775),
## num.parent.pts=5775, mode=Fixed-k, k=25, s=3e-04,
## method=TSD:vmax
## - computing values of kmax, rmax, and amax...Done
## - set of neighbors (re)named: toni|vmax|s3e-04|n5775|kmax25|rmax162.9|amax1399.3
## Finding nearest neighbors for id=toni (n=5775),
## num.parent.pts=5775, mode=Fixed-k, k=25, s=0.003,
## method=TSD:vmax
## - computing values of kmax, rmax, and amax...Done
## - set of neighbors (re)named: toni|vmax|s0.003|n5775|kmax25|rmax222.8|amax2349.7
## Finding nearest neighbors for id=toni (n=5775),
## num.parent.pts=5775, mode=Fixed-k, k=25, s=0.03, method=TSD:vmax
## - there is already a set of nearest neighbors for this set of
## parent points and value of s.
## - enough points already saved, no need to find more
## Finding nearest neighbors for id=toni (n=5775),
## num.parent.pts=5775, mode=Fixed-k, k=25, s=0.3, method=TSD:vmax
## - computing values of kmax, rmax, and amax...Done
## - set of neighbors (re)named: toni|vmax|s0.3|n5775|kmax25|rmax12947.9|amax169162.3
##
## Done. Nearest neighbor set(s) created / updated:
## toni|vmax|s3e-04|n5775|kmax25|rmax162.9|amax1399.3
## toni|vmax|s0.003|n5775|kmax25|rmax222.8|amax2349.7
## toni|vmax|s0.3|n5775|kmax25|rmax12947.9|amax169162.3
## Total time: 3 minslxy.plot.mtdr(toni.lxy, k=10)## Computing maximum theoretical distance over TSD for toni
##
|
| | 0%
|
|============= | 20%
|
|========================== | 40%
|
|======================================= | 60%
|
|==================================================== | 80%
|
|=================================================================| 100%
lxy.plot.tspan(toni.lxy, k=10)## Computing time span for toni
##
|
| | 0%
|
|============= | 20%
|
|========================== | 40%
|
|======================================= | 60%
|
|==================================================== | 80%
|
|=================================================================| 100%
These methods still leave a great deal up to the researcher, which has its pros and cons, but at least you’ve got quite a few alternatives for selecting an ideal s value for your particular purposes. We’ll stick with our 0.03 value for now, but we’ll need to create a new lhs object using this parameter:
toni.lhs.time <- lxy.lhs(toni.lxy, k=3*3:8, s=0.03)## Using nearest-neighbor selection mode: Fixed-k
## Constructing hulls and hull metrics...
## toni: 9 duplicate points were randomly displaced by 1 map unit(s)
##
## toni.pts5775.k9.s0.03.kmin0
## Found a suitable set of nearest neighbors
## Identifying the boundary points for each parent point
## Converting boundary points into polygons
## Calculating area and perimeter...Done.
## Calculating the time span of each hull...Done.
## Identifying enclosed points...Done.
##
## toni.pts5775.k12.s0.03.kmin0
## Found a suitable set of nearest neighbors
## Identifying the boundary points for each parent point
## Converting boundary points into polygons
## Calculating area and perimeter...Done.
## Calculating the time span of each hull...Done.
## Identifying enclosed points...Done.
##
## toni.pts5775.k15.s0.03.kmin0
## Found a suitable set of nearest neighbors
## Identifying the boundary points for each parent point
## Converting boundary points into polygons
## Calculating area and perimeter...Done.
## Calculating the time span of each hull...Done.
## Identifying enclosed points...Done.
##
## toni.pts5775.k18.s0.03.kmin0
## Found a suitable set of nearest neighbors
## Identifying the boundary points for each parent point
## Converting boundary points into polygons
## Calculating area and perimeter...Done.
## Calculating the time span of each hull...Done.
## Identifying enclosed points...Done.
##
## toni.pts5775.k21.s0.03.kmin0
## Found a suitable set of nearest neighbors
## Identifying the boundary points for each parent point
## Converting boundary points into polygons
## Calculating area and perimeter...Done.
## Calculating the time span of each hull...Done.
## Identifying enclosed points...Done.
##
## toni.pts5775.k24.s0.03.kmin0
## Found a suitable set of nearest neighbors
## Identifying the boundary points for each parent point
## Converting boundary points into polygons
## Calculating area and perimeter...Done.
## Calculating the time span of each hull...Done.
## Identifying enclosed points...Done.
## The following hullsets were generated:
## toni.pts5775.k9.s0.03.kmin0
## toni.pts5775.k12.s0.03.kmin0
## toni.pts5775.k15.s0.03.kmin0
## toni.pts5775.k18.s0.03.kmin0
## toni.pts5775.k21.s0.03.kmin0
## toni.pts5775.k24.s0.03.kmin0
## Total time: 26.7 secstoni.lhs.time <- lhs.iso.add(toni.lhs.time)## Merging hulls into isopleths
## toni.pts5775.k9.s0.03.kmin0
## Sorting hulls by area...Done.
## Unioning hulls
##
|
| | 0%
|
|++++++++++ | 20%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++| 100%
## toni.pts5775.k12.s0.03.kmin0
## Sorting hulls by area...Done.
## Unioning hulls
##
|
| | 0%
|
|++++++++++ | 20%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++| 100%
## toni.pts5775.k15.s0.03.kmin0
## Sorting hulls by area...Done.
## Unioning hulls
##
|
| | 0%
|
|++++++++++ | 20%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++| 100%
## 1 invalid polygon(s) removed
## toni.pts5775.k18.s0.03.kmin0
## Sorting hulls by area...Done.
## Unioning hulls
##
|
| | 0%
|
|++++++++++ | 20%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++| 100%
## toni.pts5775.k21.s0.03.kmin0
## Sorting hulls by area...Done.
## Unioning hulls
##
|
| | 0%
|
|++++++++++ | 20%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++| 100%
## toni.pts5775.k24.s0.03.kmin0
## Sorting hulls by area...Done.
## Unioning hulls
##
|
| | 0%
|
|++++++++++ | 20%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++| 100%
## Total time: 11.2 secsplot(toni.lhs.time, iso=TRUE, figs.per.page=6)
Once again, it looks like k=15 meets the eyeball-test, even when we incorporate the temporal component. We can take one last look to see what the old isopleths (without time) look like in comparison to these (with time).
toni.lhs.time.k15 <- lhs.select(toni.lhs.time, k=15)
plot(toni.lhs.k15, iso=TRUE, figs.per.page=1)
plot(toni.lhs.time.k15, iso=TRUE, figs.per.page=1)
Before we finish up this segment, let’s make sure we save our new hullset so that we can use it in the next lesson!
lhs.save(toni.lhs.time.k15)## LoCoH-hullset toni.lhs.time.k15 saved as 'toni.lhs.time.k15' to:
## /Users/ericdougherty/Desktop/GitHub/MovEco-R-Workshop/Materials/Day4/toni.n5775.s0.03.k15.iso.lhs.03.RDataOne final bit of information regarding the parameter seletion process. I am sure you have noticed that there is quite a lot of subjectivity in all of the methods we went over for selecting s and k. Though this offers a great deal of flexibility (one of the stengths of the T-LoCoH algorithm and package), it makes many researchers a bit uncomfortable. In order to circumvent some of this discomfort, several of my colleagues and I proposed an alternative method for simultaneously selecting appropriate s and k values using a grid-based search and an information criterion of our devising. Though this takes longer than the alternatives mentioned above, it does offer a fundamental theoretical underpinning to the selection process. Below is a function that wraps up this search process and outputs a data frame with values based on two different information criteria (one siliar to AIC and another akin to BIC). Let’s go through the code together, and based on what we have done above, we’ll try to figure out what this algorithm does to optimize the two parameters.
The first function train.test doesn’t actually depend on T-LoCoH at all. We are just creating a set of 100 different training/testing data sets for use in the algo function. This is where we actually create hullsets and attempt to determine the optimal parameter values. The inputs are the training/testing set we made in the previous function, a vector of the k values we want to search over, a single value for the number of different s values to test for (i.e., 40 will result in divisions of 0.025 each), the lxy object, and finally, the original dataset.
train.test <- function (data, seed) {
df <- data.frame(matrix(TRUE,nrow(data),100))
set.seed(seed)
for (i in 1:ncol(df)) {
samp <- sample(seq(1,nrow(data),1), round(0.002222*(nrow(data))))
for (j in 1:length(samp)) {
for (k in 1:nrow(df)) {
if (k == samp[j]) {
df[k,i] <- FALSE
}
}
}
}
count <- 0
for (i in 1:ncol(df)) {
new <- table(df[,i])[1]
count <- count + new
}
for (i in 3:nrow(df)) {
for (j in 1:ncol(df)) {
if (df[i-1,j] == FALSE && df[i-2,j] != FALSE) {
for (k in 0:98) {
df[i+k,j] <- FALSE
}
}
}
}
df <- df[1:nrow(data),]
return(df)
}
algo <- function(train.test, k.vals, num.s.vals, data.lxy, data) {
# Determine number of test points across all 100 train/test splits for BIC calculation
count <- 0
for (i in 1:ncol(train.test)) {
new <- table(train.test[,i])[1]
count <- count + new
}
trace <- data.frame(matrix(0,length(k.vals),7))
for (k in 1:length(k.vals)) {
for (z in 1:(num.s.vals + 1)) {
current.k_val <- k.vals[k]
current.s_val <- (z-1) * (1/num.s.vals)
# Calculate the nearest neighbors and create lhs object for full dataset
full.lxy <- lxy.nn.add(data.lxy, k=current.k_val, s=current.s_val, status=F)
full.lhs <- lxy.lhs(full.lxy, k=current.k_val, s=current.s_val, status=F)
coords <- full.lhs[[1]]$pts@coords
# Create list for the negative log likelihood values for each test/train split in df
likelihood <- list()
for (j in 1:ncol(train.test)) {
# Create a one-column data frame from df
df1 <- df[1:nrow(data),j]
# Create selection of hulls based on Boolean
hulls.sel.idx <- which(df1)
full.hulls <- hulls(full.lhs)
selected.hulls <- full.hulls[[1]] [ full.hulls[[1]]@data$pts.idx %in% hulls.sel.idx , ]
# Determine the number of points in training dataset
df.temp <- data.frame(as.numeric(train.test[,j]))
colnames(df.temp) <- "subset"
total.pts <- sum(df.temp)
# Find middle points of testing data and define as -1
for (i in 2:nrow(df.temp)) {
if (df.temp[i,1] == 0 && df.temp[i-1,1] != 0 && df.temp[i-1,1] != -1) {
df.temp[i+50,1] <- -1
}
}
df.temp <- df.temp[1:nrow(data),]
df.temp <- cbind(coords, df.temp)
# Extract the middle points for testing and record associated coordinates
test.pts <- data.frame()
q = 1
for (i in 1:nrow(df.temp)) {
if (df.temp[i,3] == -1) {
test.pts[q,1] <- df.temp[i,1]
test.pts[q,2] <- df.temp[i,2]
q = q + 1
}
}
colnames(test.pts) <- c("x", "y")
test.pts <- SpatialPoints(test.pts[ , c("x","y")], proj4string=CRS("+proj=utm +south +zone=35 +ellps=WGS84"))
poly <- SpatialPolygons(selected.hulls@polygons, proj4string = CRS('+proj=utm +south +zone=35 +ellps=WGS84'))
# Determine the number of hulls under each test point,
overlay <- data.frame(matrix(0,length(test.pts@coords[,1]),1))
for (i in 1:length(test.pts@coords[,1])) {
overlay.list <- over(test.pts[i], poly, returnList=TRUE)
overlay[i,1] <- length(overlay.list[[1]])
}
hull.mean <- mean(full.lhs[[1]]$hulls@data$area)
# Calculate likelihood
for (i in 1:nrow(overlay)) {
overlay[i,2] <- overlay[i,1]/length(selected.hulls[[1]])
overlay[i,3] <- -log(overlay[i,2])
overlay[i,4] <- log(overlay[i,2],exp(1))
if (overlay[i,1] == 0) {
overlay[i,3] <- -log((1/nrow(data))/100)
overlay[i,4] <- log((1/nrow(data))/100, exp(1))
overlay[i,4] <- log((1/nrow(data))/100, exp(1))
}
}
colnames(overlay) <- c("over", "prob", "loglike", "ln.like")
# Add values to likelihood list
likelihood[[j]] <- as.list(overlay)
}
log.like <- data.frame(matrix(0,length(likelihood),4))
for (i in 1:length(likelihood)) {
log.like[i,1] <- sum(likelihood[[i]]$loglike, na.rm=TRUE)
log.like[i,2] <- sum(likelihood[[i]]$ln.like, na.rm=TRUE)
log.like[i,3] <- -2*(log.like[i,2]) + 2*k
log.like[i,4] <- -2*(log.like[i,2]) + k*log(as.numeric(count),exp(1))
}
# Assign mean of the means of negative log likelihoods as the current posterior probability for Metropolis-Hastings
new.postLike <- sum(log.like[,1])
new.lnLike <- sum(log.like[,2])
new.AIC <- sum(log.like[,3])
new.BIC <- sum(log.like[,4])
#cat("iteration:", counter, "chain:", current.k_val, current.s_val, "likelihood:", new.postLike, "\n")
trace[k,1] <- current.s_val
trace[k,2] <- current.k_val
trace[k,3] <- hull.mean
trace[k,4] <- new.postLike
trace[k,5] <- new.lnLike
trace[k,6] <- new.AIC
trace[k,7] <- new.BIC
} # End of s loop
} # End of k loop
return(trace)
}Just to show what the result would look like without going through the multi-hour computation, we can call in an existing trace output:
trace <- read.csv('Example_Trace.csv')
#install.packages('lattice')
library(lattice)## Warning: package 'lattice' was built under R version 3.3.2my.palette <- c("#FFF5F0", "#FEE0D2", "#FEE0D2", "#FCBBA1", "#FCBBA1", "#FC9272", "#FC9272", "#FB6A4A", "#FB6A4A", "#EF3B2C", "#EF3B2C", "#CB181D", "#CB181D", "#A50F15", "#A50F15", "#67000D")
levelplot(-trace$X7 ~ trace$X2 * trace$X3, xlab="k value", ylab="s value", zlab="IC", screen = list(z = -30, x=-60), regions=TRUE, cuts=15, col.regions=my.palette)
wireframe(-trace$X7 ~ trace$X2 * trace$X3, xlab="k value", ylab="s value", zlab="IC", screen = list(z = -30, x=-60), drape=TRUE, cuts=15, col.regions=my.palette)